%Barmak's method for calculate Zernike radial polynomials
function RM = Barmak(j_max,n_max,lent,r)
%RM(1,1) = 1; %R_00
%独立计算n=m的时候
RM = zeros(j_max-1,lent);
RM(1:2,:) = [r;r];

for n = 2:n_max %calculate cases when n=m
    j = (n^2+3*n)/2;
    j_pre = ((n-1)^2+3*(n-1))/2;
    RM(j-1:j,:) = [r;r].*RM(j_pre-1:j_pre,:);
end

RM(3,:) = 2*r.*(RM(1,:)) - 1; 
%RM(3,:) = 2*r.*(RM(1,:)) - 0;
for n = 3:n_max %n
    for m = 0:n-2  %m
        if rem((n-m),2) == 0
            index_r = [n-2,m];
            index_h = [n-1,m+1];
            index_b = [n-1,m-1,m+1];
            if m == 0
                j = n*(n+1)/2+1;
                j_r = index_r(1)*(index_r(1)+1)/2+1;
                j_h = index_h(1)*(index_h(1)+1)/2+index_h(2);
                RM(j-1,:) = 2*r.*RM(j_h-1,:) - RM(j_r-1,:);
            elseif m == 1    
                j = n*(n+1)/2+m;
                j_r = index_r(1)*(index_r(1)+1)/2+index_r(2);
                j_b1 = index_b(1)*(index_b(1)+1)/2+1;
                j_b2 = index_b(1)*(index_b(1)+1)/2+index_b(3);
                RM(j-1:j,:) = repmat(r.*(RM(j_b1-1,:)+RM(j_b2-1,:))...
                    - RM(j_r-1,:),2,1);
            else
                j = n*(n+1)/2+m;
                j_r = index_r(1)*(index_r(1)+1)/2+index_r(2);
                j_b1 = index_b(1)*(index_b(1)+1)/2+index_b(2);
                j_b2 = index_b(1)*(index_b(1)+1)/2+index_b(3);
                RM(j-1:j,:) = [r;r].*(RM(j_b1-1:j_b1,:)+RM(j_b2-1:j_b2,:))...
                    - RM(j_r-1:j_r,:);
            end
        end
    end
    
end
end